Chadha et al. (2019) set up a completely randomized experiment studying the effects of different soil water holding capacities (WHCs) on various characteristics of the agronomic weed Lactuca serriola. Plants were grown individually in pots allocated to one of four WHCs (100%, 75%, 50% and 25%). There were seven plants for each WHC, although we just use a subset of the data from the 100% treatment group so there was no between-subjects factor. The number of leaves on each of the seven plants in 100% WHC soil was recorded weekly for nine weeks (counts at the start of the experiment, week 0, were omitted). Time was the within-subjects (repeated measures) fixed factor and individual plants were the random subjects. With nine weeks and reasonably linear trends through time for each plant (Figure 12.2), it made sense to treat time as a continuous covariate for analysis. Although week 0 was not included in the analysis, we did not centre time for analysis, so intercepts represent the number of leaves for week 0.

Prickly lettuce or Milk thistle. Mick Keough, CC SA-BY 4.0

Chadha, A., Florentine, S., Chauhan, B. S., Long, B. & Jayasundera, M. (2019). Influence of soil moisture regimes on growth, photosynthetic capacity, leaf biochemistry and reproductive capabilities of the invasive agronomic weed; Lactuca serriola. PLoS One, 14, e0218191.

Link to paper

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans, Rmisc, MuMIn)

Import chadha data file

chadha <- read.csv("../data/chadha.csv")
chadha

set contrasts from afex

set_sum_contrasts()
## setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))

select 100%WHC

chadha1 <- subset(chadha, treat=="100WHC")

make individual a factor

chadha1$plant <- factor(chadha1$plant)

plot slopes - very consistent

xyplot(noleaves~week|plant, type=c("p","r"), auto.key=T, chadha1)

Fit OLS “ancova” model

chadha1.aov <- aov(noleaves~week*plant, data=chadha1)
plot(chadha1.aov)

summary(chadha1.aov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## week         1   4999    4999 310.550  < 2e-16 ***
## plant        6    710     118   7.352 1.22e-05 ***
## week:plant   6    110      18   1.143    0.352    
## Residuals   49    789      16                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Get Type III SS)

Anova(lm(chadha1.aov), type='III')

Get correct F-ratio and P value for week (tested against week by plant)

4999.050/18.400
## [1] 271.6875
1-pf(4999.050/18.400, 1, 6, lower.tail = TRUE, log.p = FALSE)
## [1] 3.178231e-06

no GG and HF adjustments as week is continuous with 1 df

Get variance components

chadha1.vca <- anovaMM(noleaves~week+(plant)+(week*plant), chadha1)
chadha1.vca
## 
## 
## ANOVA-Type Estimation of Mixed Model:
## --------------------------------------
## 
##  [Fixed Effects]
## 
##      int     week 
## 16.40079  3.45000 
## 
## 
##  [Variance Components]
## 
##   Name       DF        SS         MS         VC        %Total    SD      
## 1 total      15.101789                       26.537662 100       5.151472
## 2 plant      6         710.095238 118.349206 10.401893 39.196719 3.225197
## 3 week:plant 6         110.4      18.4       0.038377  0.144613  0.1959  
## 4 error      49        788.772222 16.097392  16.097392 60.658668 4.012156
##   CV[%]    
## 1 15.308619
## 2 9.584311 
## 3 0.582156 
## 4 11.922915
## 
## Mean: 33.65079 (N = 63) 
## 
## Experimental Design: balanced  |  Method: ANOVA
VCAinference(chadha1.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)
## 
## 
## 
## Inference from Mixed Model Fit
## ------------------------------
## 
## > VCA Result:
## -------------
## 
##  [Fixed Effects]
## 
##     int    week 
## 16.4008  3.4500 
## 
## 
##  [Variance Components]
## 
##   Name       DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
## 1 total      15.1018                   26.5377 100     5.1515 15.3086 92.6777
## 2 plant      6       710.0952 118.3492 10.4019 39.1967 3.2252 9.5843  77.6065
## 3 week:plant 6       110.4    18.4     0.0384  0.1446  0.1959 0.5822  0.0343 
## 4 error      49      788.7722 16.0974  16.0974 60.6587 4.0122 11.9229 10.5766
## 
## Mean: 33.6508 (N = 63) 
## 
## Experimental Design: balanced  |  Method: ANOVA
## 
## 
## > VC:
## -----
##            Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
## total       26.5377 14.5063 63.3391       15.9486       54.6613
## plant       10.4019 -6.8643 27.6681       -4.0884       24.8922
## week:plant   0.0384 -0.3245  0.4013       -0.2662        0.3429
## error       16.0974 11.2325 24.9968       11.8901       23.2468
## 
## > SD:
## -----
##            Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
## total        5.1515  3.8087 7.9586        3.9936        7.3933
## plant        3.2252 -2.6200 5.2600       -2.0220        4.9892
## week:plant   0.1959 -0.5697 0.6335       -0.5159        0.5856
## error        4.0122  3.3515 4.9997        3.4482        4.8215
## 
## > CV[%]:
## --------
##            Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
## total       15.3086 11.3183 23.6505       11.8677       21.9707
## plant        9.5843 -7.7858 15.6313       -6.0087       14.8264
## week:plant   0.5822 -1.6929  1.8825       -1.5332        1.7403
## error       11.9229  9.9596 14.8575       10.2470       14.3280
## 
## 
## 95% Confidence Level  
## SAS PROC MIXED method used for computing CIs

Drop interaction

chadha2.aov <- aov(noleaves~week+plant, data=chadha1)
Anova(lm(chadha2.aov), type='III')

get new variance components

chadha2.vca <- anovaMM(noleaves~week+(plant), chadha1)
chadha2.vca
## 
## 
## ANOVA-Type Estimation of Mixed Model:
## --------------------------------------
## 
##  [Fixed Effects]
## 
##      int     week 
## 16.40079  3.45000 
## 
## 
##  [Variance Components]
## 
##   Name  DF        SS         MS         VC        %Total    SD       CV[%]    
## 1 total 23.462944                       27.681988 100       5.261368 15.635196
## 2 plant 6         710.095238 118.349206 11.333402 40.941432 3.366512 10.004257
## 3 error 55        899.172222 16.348586  16.348586 59.058568 4.043338 12.015581
## 
## Mean: 33.65079 (N = 63) 
## 
## Experimental Design: balanced  |  Method: ANOVA
VCAinference(chadha2.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)
## 
## 
## 
## Inference from Mixed Model Fit
## ------------------------------
## 
## > VCA Result:
## -------------
## 
##  [Fixed Effects]
## 
##     int    week 
## 16.4008  3.4500 
## 
## 
##  [Variance Components]
## 
##   Name  DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
## 1 total 23.4629                   27.682  100     5.2614 15.6352 65.3194
## 2 plant 6       710.0952 118.3492 11.3334 40.9414 3.3665 10.0043 57.76  
## 3 error 55      899.1722 16.3486  16.3486 59.0586 4.0433 12.0156 9.7191 
## 
## Mean: 33.6508 (N = 63) 
## 
## Experimental Design: balanced  |  Method: ANOVA
## 
## 
## > VC:
## -----
##       Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
## total  27.6820 16.7947 54.0455       18.1687       48.3233
## plant  11.3334 -3.5623 26.2291       -1.1675       23.8343
## error  16.3486 11.6201 24.7038       12.2651       23.0805
## 
## > SD:
## -----
##       Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
## total   5.2614  4.0981 7.3516        4.2625        6.9515
## plant   3.3665 -1.8874 5.1214       -1.0805        4.8820
## error   4.0433  3.4088 4.9703        3.5022        4.8042
## 
## > CV[%]:
## --------
##       Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
## total  15.6352 12.1784 21.8466       12.6668       20.6578
## plant  10.0043 -5.6088 15.2194       -3.2109       14.5079
## error  12.0156 10.1300 14.7702       10.4073       14.2767
## 
## 
## 95% Confidence Level  
## SAS PROC MIXED method used for computing CIs

Now mixed effects modelcomparing random slopes and random intercepts using ML

chadha.lmer1 <- lmer(noleaves~week + (week|plant), REML=FALSE, chadha1)
## boundary (singular) fit: see help('isSingular')
 # singular fit so set correlation between random effects to zero
chadha.lmer1 <- lmer(noleaves~week + (week||plant), REML=FALSE, chadha1)
chadha.lmer2 <- lmer(noleaves~week + (1|plant), REML=FALSE, chadha1)

Compare models

anova(chadha.lmer1, chadha.lmer2)
AICc(chadha.lmer1, chadha.lmer2)

Focus on random intercept model - refit with REML

chadha.lmer3 <- lmer(noleaves~week + (1|plant), REML=TRUE, chadha1)
summary(chadha.lmer3, ddf="Kenward-Roger")
## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: noleaves ~ week + (1 | plant)
##    Data: chadha1
## 
## REML criterion at convergence: 365.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9226 -0.5086  0.0859  0.3199  3.9265 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plant    (Intercept) 11.33    3.367   
##  Residual             16.35    4.043   
## Number of obs: 63, groups:  plant, 7
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  16.4008     1.6887 13.4331   9.712 1.91e-07 ***
## week          3.4500     0.1973 55.0000  17.487  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## week -0.584
anova(chadha.lmer3, type=3, ddf="Kenward-Roger")

Get variance components

chadha.ci3 <- confint.merMod(chadha.lmer3, oldNames=FALSE)
## Computing profile confidence intervals ...
chadha.vc3 <- (chadha.ci3)^2
print(chadha.vc3)
##                           2.5 %    97.5 %
## sd_(Intercept)|plant   2.765720  38.24037
## sigma                 11.327647  23.82822
## (Intercept)          170.118317 390.40409
## week                   9.364274  14.74475

Random slopes output for comparison - refit using REML

chadha.lmer4 <- lmer(noleaves~week + (week||plant), REML=TRUE, chadha1)
summary(chadha.lmer4, ddf="Kenward-Roger")
## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: noleaves ~ week + (week || plant)
##    Data: chadha1
## 
## REML criterion at convergence: 364.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0166 -0.5212  0.0766  0.3747  3.9247 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plant    (Intercept)  6.797   2.6071  
##  plant.1  week         0.141   0.3754  
##  Residual             15.644   3.9553  
## Number of obs: 63, groups:  plant, 7
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  16.4008     1.4665 10.2815   11.18 4.41e-07 ***
## week          3.4500     0.2396 10.2815   14.40 3.79e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## week -0.530
anova(chadha.lmer4, type=3, ddf="Kenward-Roger")

Check whether AR(1) covariance structure would improve the fit using nlme

chadha.lme1 <- lme(noleaves~week, random=~1|plant, method="ML", chadha1)
chadha.lme2 <- lme(noleaves~week, random=~1|plant, method="ML", correlation=corAR1(form=~1|plant), chadha1)
anova(chadha.lme1, chadha.lme2)
AICc(chadha.lme1, chadha.lme2)

no improvement with AR(1)

chadha <- lm(noleaves~week, chadha1)
summary(chadha)
## 
## Call:
## lm(formula = noleaves ~ week, data = chadha1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1008 -3.1508  0.2492  2.9492 14.5492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  16.4008     1.4103   11.63   <2e-16 ***
## week          3.4500     0.2506   13.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.136 on 61 degrees of freedom
## Multiple R-squared:  0.7565, Adjusted R-squared:  0.7525 
## F-statistic: 189.5 on 1 and 61 DF,  p-value: < 2.2e-16